GLMM example 3

Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(brms)
library(tidybayes)
library(broom.mixed) # for tidying MCMC outputs
library(patchwork) # for multiple plots
library(standist) # for visualizing distributions
library(rstanarm)
library(ggeffects)
library(bayesplot)
library(rstan)
library(DHARMa)
source("helperFunctions.R")

Scenario

Starlings

Format of starling_full.RSV data files

SITUATION	MONTH	MASS	BIRD
tree	Nov	78	tree1
..	..	..	..
nest-box	Nov	78	nest-box1
..	..	..	..
inside	Nov	79	inside1
..	..	..	..
other	Nov	77	other1
..	..	..	..
tree	Jan	85	tree1
..	..	..	..

SITUATION	Categorical listing of roosting situations (tree, nest-box, inside or other)
MONTH	Categorical listing of the month of sampling.
MASS	Mass (g) of starlings.
BIRD	Categorical listing of individual bird repeatedly sampled.

Read in the data

starling <- read_csv("../public/data/starling_full.csv", trim_ws = TRUE)
starling %>% glimpse()

## Rows: 80
## Columns: 5
## $ MONTH      <chr> "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "No…
## $ SITUATION  <chr> "tree", "tree", "tree", "tree", "tree", "tree", "tree", "tr…
## $ subjectnum <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1…
## $ BIRD       <chr> "tree1", "tree2", "tree3", "tree4", "tree5", "tree6", "tree…
## $ MASS       <dbl> 78, 88, 87, 88, 83, 82, 81, 80, 80, 89, 78, 78, 85, 81, 78,…

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i = \beta_0 + \boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \boldsymbol{\gamma} = \gamma_0\\ \beta_0 \sim{} \mathcal{N}(0, 100)\\ \beta \sim{} \mathcal{N}(0, 10)\\ \gamma_0 \sim{} \mathcal{N}(0, \sigma_1^2)\\ \sigma \sim{} \mathcal{cauchy}(0, 2)\\ \sigma_1 \sim{} \mathcal{cauchy}(0, 2)\\ \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of roosting situation and month on starling mass. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual birds.

starling <- starling %>%
  mutate(
    BIRD = factor(BIRD),
    SITUATION = factor(SITUATION),
    MONTH = factor(MONTH, levels = c("Nov", "Jan"))
  )

ggplot(starling, aes(y = MASS, x = MONTH)) +
  geom_boxplot() +
  facet_grid(~SITUATION)

## Better still
ggplot(starling, aes(y = MASS, x = MONTH, group = BIRD)) +
  geom_point() +
  geom_line() +
  facet_grid(~SITUATION)

Conclusions:

• it is clear that the Nov mass of each bird is different - so random intercepts
• the degree to which they change between Nov and Dec is also relatively different per bird - perhaps random intercept/random slope

Fit the model

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

starling.rstanarm <- stan_glmer(MASS ~ MONTH * SITUATION + (1 | BIRD),
  data = starling,
  family = gaussian(),
  iter = 5000,
  warmup = 2000,
  chains = 3,
  thin = 5,
  refresh = 0
)

starling.rstanarm %>% prior_summary()

## Priors for model '.' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 84, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 84, scale = 17)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [33.25,38.40,38.40,...])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.15)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

This tells us: - for the intercept (when the family is Gaussian), a normal prior with a mean of 84 and a standard deviation of 17 is used. The 2.5 is used for scaling all parameter standard deviations. The value of 84 is based on the mean of the response (mean(starling$MASS)) and the scaled standard deviation of 17 is based on multiplying the scaling factor by the standard deviation of the response.

2.5 * sd(starling$MASS)

## [1] 16.73225

• for the coefficients (in this case, the effects of MONTH, SEASON and their interactions), the default prior is a normal prior centred around 0 with a standard deviation of 2.5. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the standard deviation of the numerical dummy variables for the predictor (then rounded).

2.5 * sd(starling$MASS) / apply(model.matrix(~ MONTH * SITUATION, starling)[, -1], 2, sd)

##                   MONTHJan          SITUATIONnest-box 
##                   33.25470                   38.39922 
##             SITUATIONother              SITUATIONtree 
##                   38.39922                   38.39922 
## MONTHJan:SITUATIONnest-box    MONTHJan:SITUATIONother 
##                   50.27638                   50.27638 
##     MONTHJan:SITUATIONtree 
##                   50.27638

• the auxiliary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be \(\sigma\), for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in rstanarm, this is a exponential with a rate of 1 which is then adjusted by division with the standard deviation of the response.

1 / sd(starling$MASS)

## [1] 0.149412

Assessing priors

Lets now run with priors only so that we can explore the range of values they allow in the posteriors.

starling.rstanarm1 <- update(starling.rstanarm, prior_PD = TRUE)

starling.rstanarm1 %>%
  ggpredict(~ MONTH * SITUATION) %>%
  plot(add.data = TRUE)

Conclusions:

• we see that the range of predictions is fairly wide and the predicted means could range from a small value to a relatively large value.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

• \(\beta_0\): normal centred at 84 with a standard deviation of 17
  • mean of 84: since mean(starling$MASS)
  • sd of 17: since 2.5*sd(starling$MASS)
• \(\beta_1\): normal centred at 0 with a standard deviation of 33
  • sd of 33: since 2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)
• \(\beta_{2-4}\): normal centred at 0 with a standard deviation of 39
  • sd of 39: since 2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)
• \(\beta_{5-7}\): normal centred at 0 with a standard deviation of 50
  • sd of 50: since 2.5*sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[,-1], 2, sd)
• \(\sigma\): exponential with rate of 0.15
  • since 1 / sd(starling$MASS)
• \(\Sigma\): decov with:
  • regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

I will also overlay the raw data for comparison.

starling.rstanarm2 <- stan_glmer(MASS ~ MONTH * SITUATION + (1 | BIRD),
  data = starling,
  family = gaussian(),
  prior_intercept = normal(84, 17, autoscale = FALSE),
  prior = normal(0, c(33, 39, 39, 39, 50, 50, 50), autoscale = FALSE),
  prior_aux = rstanarm::exponential(0.15, autoscale = FALSE),
  prior_covariance = decov(1, 1, 1, 1),
  prior_PD = TRUE,
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

starling.rstanarm2 %>%
  ggpredict(~ SITUATION * MONTH) %>%
  plot(add.data = TRUE)

Now lets refit, conditioning on the data.

starling.rstanarm3 <- update(starling.rstanarm2, prior_PD = FALSE)

Plotting prior and posterior

posterior_vs_prior(starling.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggemmeans(starling.rstanarm3, ~ SITUATION * MONTH) %>% plot(add.data = TRUE)

ggpredict(starling.rstanarm3, ~ SITUATION * MONTH) %>% plot(add.data = TRUE)

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

starling.form <- bf(MASS ~ MONTH*SITUATION+(1|BIRD),
                   family = gaussian() 
                   )
options(width=150)
starling.form %>% get_prior(data = starling)

##                  prior     class                       coef group resp dpar nlpar bound       source
##                 (flat)         b                                                             default
##                 (flat)         b                   MONTHJan                             (vectorized)
##                 (flat)         b MONTHJan:SITUATIONnestMbox                             (vectorized)
##                 (flat)         b    MONTHJan:SITUATIONother                             (vectorized)
##                 (flat)         b     MONTHJan:SITUATIONtree                             (vectorized)
##                 (flat)         b          SITUATIONnestMbox                             (vectorized)
##                 (flat)         b             SITUATIONother                             (vectorized)
##                 (flat)         b              SITUATIONtree                             (vectorized)
##  student_t(3, 84, 5.9) Intercept                                                             default
##   student_t(3, 0, 5.9)        sd                                                             default
##   student_t(3, 0, 5.9)        sd                             BIRD                       (vectorized)
##   student_t(3, 0, 5.9)        sd                  Intercept  BIRD                       (vectorized)
##   student_t(3, 0, 5.9)     sigma                                                             default

options(width=80)

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

• \(\beta_0\): normal centred at 84 with a standard deviation of 10
• \(\beta_1\): normal centred at 0 with a standard deviation of 10
• \(\sigma\): gamma with shape parameters 2 and 0.5
• \(\sigma_j\): half-cauchy with parameters 0 and 5.
• \(\Sigma\): decov with:
  • regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

• half-t: as the degrees of freedom approach infinity, this will approach a half-normal
• half-cauchy: this is essentially a half-t with 1 degree of freedom
• exponential

I will also overlay the raw data for comparison.

starling %>%
  group_by(SITUATION, MONTH) %>%
  summarise(
    median(MASS),
    mad(MASS)
  )

standist::visualize("normal(84,20)", xlim = c(0, 200))

standist::visualize("student_t(3, 0, 5.9)",
  "gamma(2,0.5)",
  "cauchy(0,1)",
  xlim = c(-10, 25)
)

priors <- prior(normal(84, 5), class = "Intercept") +
  prior(normal(0, 5), class = "b") +
  prior(gamma(2, 0.5), class = "sigma") +
  prior(cauchy(0, 1), class = "sd")
starling.form <- bf(MASS ~ MONTH * SITUATION + (1 | BIRD),
  family = gaussian()
)
starling.brm2 <- brm(starling.form,
  data = starling,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

starling.form <- bf(MASS ~ MONTH * SITUATION + (MONTH | BIRD),
  family = gaussian()
)
starling.brm3 <- brm(starling.form,
  data = starling,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0,
  control = list(adapt_delta = 0.99)
)

(l.1 <- starling.brm2 %>% loo())

## 
## Computed from 2400 by 80 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -233.0  5.1
## p_loo        10.8  1.3
## looic       466.1 10.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     74    92.5%   1088      
##  (0.5, 0.7]   (ok)        6     7.5%   642       
##    (0.7, 1]   (bad)       0     0.0%   <NA>      
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

(l.2 <- starling.brm3 %>% loo())

## 
## Computed from 2400 by 80 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -232.6  5.1
## p_loo        15.4  1.8
## looic       465.2 10.3
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     54    67.5%   825       
##  (0.5, 0.7]   (ok)       22    27.5%   196       
##    (0.7, 1]   (bad)       4     5.0%   1021      
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.

loo_compare(l.1, l.2)

##   elpd_diff se_diff
## .  0.0       0.0   
## . -0.4       0.5

starling.brm3 %>% get_variables()

##   [1] "b_Intercept"                   "b_MONTHJan"                   
##   [3] "b_SITUATIONnestMbox"           "b_SITUATIONother"             
##   [5] "b_SITUATIONtree"               "b_MONTHJan:SITUATIONnestMbox" 
##   [7] "b_MONTHJan:SITUATIONother"     "b_MONTHJan:SITUATIONtree"     
##   [9] "sd_BIRD__Intercept"            "sd_BIRD__MONTHJan"            
##  [11] "cor_BIRD__Intercept__MONTHJan" "sigma"                        
##  [13] "r_BIRD[inside1,Intercept]"     "r_BIRD[inside10,Intercept]"   
##  [15] "r_BIRD[inside2,Intercept]"     "r_BIRD[inside3,Intercept]"    
##  [17] "r_BIRD[inside4,Intercept]"     "r_BIRD[inside5,Intercept]"    
##  [19] "r_BIRD[inside6,Intercept]"     "r_BIRD[inside7,Intercept]"    
##  [21] "r_BIRD[inside8,Intercept]"     "r_BIRD[inside9,Intercept]"    
##  [23] "r_BIRD[nest-box1,Intercept]"   "r_BIRD[nest-box10,Intercept]" 
##  [25] "r_BIRD[nest-box2,Intercept]"   "r_BIRD[nest-box3,Intercept]"  
##  [27] "r_BIRD[nest-box4,Intercept]"   "r_BIRD[nest-box5,Intercept]"  
##  [29] "r_BIRD[nest-box6,Intercept]"   "r_BIRD[nest-box7,Intercept]"  
##  [31] "r_BIRD[nest-box8,Intercept]"   "r_BIRD[nest-box9,Intercept]"  
##  [33] "r_BIRD[other1,Intercept]"      "r_BIRD[other10,Intercept]"    
##  [35] "r_BIRD[other2,Intercept]"      "r_BIRD[other3,Intercept]"     
##  [37] "r_BIRD[other4,Intercept]"      "r_BIRD[other5,Intercept]"     
##  [39] "r_BIRD[other6,Intercept]"      "r_BIRD[other7,Intercept]"     
##  [41] "r_BIRD[other8,Intercept]"      "r_BIRD[other9,Intercept]"     
##  [43] "r_BIRD[tree1,Intercept]"       "r_BIRD[tree10,Intercept]"     
##  [45] "r_BIRD[tree2,Intercept]"       "r_BIRD[tree3,Intercept]"      
##  [47] "r_BIRD[tree4,Intercept]"       "r_BIRD[tree5,Intercept]"      
##  [49] "r_BIRD[tree6,Intercept]"       "r_BIRD[tree7,Intercept]"      
##  [51] "r_BIRD[tree8,Intercept]"       "r_BIRD[tree9,Intercept]"      
##  [53] "r_BIRD[inside1,MONTHJan]"      "r_BIRD[inside10,MONTHJan]"    
##  [55] "r_BIRD[inside2,MONTHJan]"      "r_BIRD[inside3,MONTHJan]"     
##  [57] "r_BIRD[inside4,MONTHJan]"      "r_BIRD[inside5,MONTHJan]"     
##  [59] "r_BIRD[inside6,MONTHJan]"      "r_BIRD[inside7,MONTHJan]"     
##  [61] "r_BIRD[inside8,MONTHJan]"      "r_BIRD[inside9,MONTHJan]"     
##  [63] "r_BIRD[nest-box1,MONTHJan]"    "r_BIRD[nest-box10,MONTHJan]"  
##  [65] "r_BIRD[nest-box2,MONTHJan]"    "r_BIRD[nest-box3,MONTHJan]"   
##  [67] "r_BIRD[nest-box4,MONTHJan]"    "r_BIRD[nest-box5,MONTHJan]"   
##  [69] "r_BIRD[nest-box6,MONTHJan]"    "r_BIRD[nest-box7,MONTHJan]"   
##  [71] "r_BIRD[nest-box8,MONTHJan]"    "r_BIRD[nest-box9,MONTHJan]"   
##  [73] "r_BIRD[other1,MONTHJan]"       "r_BIRD[other10,MONTHJan]"     
##  [75] "r_BIRD[other2,MONTHJan]"       "r_BIRD[other3,MONTHJan]"      
##  [77] "r_BIRD[other4,MONTHJan]"       "r_BIRD[other5,MONTHJan]"      
##  [79] "r_BIRD[other6,MONTHJan]"       "r_BIRD[other7,MONTHJan]"      
##  [81] "r_BIRD[other8,MONTHJan]"       "r_BIRD[other9,MONTHJan]"      
##  [83] "r_BIRD[tree1,MONTHJan]"        "r_BIRD[tree10,MONTHJan]"      
##  [85] "r_BIRD[tree2,MONTHJan]"        "r_BIRD[tree3,MONTHJan]"       
##  [87] "r_BIRD[tree4,MONTHJan]"        "r_BIRD[tree5,MONTHJan]"       
##  [89] "r_BIRD[tree6,MONTHJan]"        "r_BIRD[tree7,MONTHJan]"       
##  [91] "r_BIRD[tree8,MONTHJan]"        "r_BIRD[tree9,MONTHJan]"       
##  [93] "prior_Intercept"               "prior_b"                      
##  [95] "prior_sigma"                   "prior_sd_BIRD"                
##  [97] "prior_cor_BIRD"                "lp__"                         
##  [99] "accept_stat__"                 "stepsize__"                   
## [101] "treedepth__"                   "n_leapfrog__"                 
## [103] "divergent__"                   "energy__"

starling.brm3 %>%
  hypothesis("MONTHJan=0") %>%
  plot()

starling.brm3 %>%
  posterior_samples() %>%
  dplyr::select(-`lp__`) %>%
  pivot_longer(everything(), names_to = "key") %>%
  filter(!str_detect(key, "^r")) %>%
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    ## Class = ifelse(str_detect(key, 'Intercept'),  'Intercept',
    ##         ifelse(str_detect(key, 'b'),  'b', 'sigma')),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_SITUATION.*|b_MONTH.*|prior_b") ~ "TREATMENT",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) %>%
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
  facet_wrap(~Class, scales = "free")

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

pars <- starling.brm3 %>% get_variables()
pars <- pars %>%
  str_extract("^b.Intercept|^b_SITUTATION.*|^b_MONTH.*|[sS]igma|^sd.*") %>%
  na.omit()
pars

## [1] "b_Intercept"                  "b_MONTHJan"                  
## [3] "b_MONTHJan:SITUATIONnestMbox" "b_MONTHJan:SITUATIONother"   
## [5] "b_MONTHJan:SITUATIONtree"     "sd_BIRD__Intercept"          
## [7] "sd_BIRD__MONTHJan"            "sigma"                       
## [9] "sigma"                       
## attr(,"na.action")
##  [1]   3   4   5  11  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27
## [20]  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46
## [39]  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65
## [58]  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84
## [77]  85  86  87  88  89  90  91  92  93  94  96  97  98  99 100 101 102 103 104
## attr(,"class")
## [1] "omit"

starling.brm3 %>% mcmc_plot(type = "trace", pars = pars)

The chains appear well mixed and very similar

• acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

starling.brm3 %>% mcmc_plot(type = "acf_bar", pars = pars)

There is no evidence of auto-correlation in the MCMC samples

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

starling.brm3 %>% mcmc_plot(type = "rhat_hist")

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

starling.brm2 %>% mcmc_plot(type = "neff_hist")

Ratios all very high.

More diagnostics

starling.brm3 %>% mcmc_plot(type = "combo", pars = pars)

starling.brm3 %>% mcmc_plot(type = "violin", pars = pars)

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

starling.brm3 %>% get_variables()

##   [1] "b_Intercept"                   "b_MONTHJan"                   
##   [3] "b_SITUATIONnestMbox"           "b_SITUATIONother"             
##   [5] "b_SITUATIONtree"               "b_MONTHJan:SITUATIONnestMbox" 
##   [7] "b_MONTHJan:SITUATIONother"     "b_MONTHJan:SITUATIONtree"     
##   [9] "sd_BIRD__Intercept"            "sd_BIRD__MONTHJan"            
##  [11] "cor_BIRD__Intercept__MONTHJan" "sigma"                        
##  [13] "r_BIRD[inside1,Intercept]"     "r_BIRD[inside10,Intercept]"   
##  [15] "r_BIRD[inside2,Intercept]"     "r_BIRD[inside3,Intercept]"    
##  [17] "r_BIRD[inside4,Intercept]"     "r_BIRD[inside5,Intercept]"    
##  [19] "r_BIRD[inside6,Intercept]"     "r_BIRD[inside7,Intercept]"    
##  [21] "r_BIRD[inside8,Intercept]"     "r_BIRD[inside9,Intercept]"    
##  [23] "r_BIRD[nest-box1,Intercept]"   "r_BIRD[nest-box10,Intercept]" 
##  [25] "r_BIRD[nest-box2,Intercept]"   "r_BIRD[nest-box3,Intercept]"  
##  [27] "r_BIRD[nest-box4,Intercept]"   "r_BIRD[nest-box5,Intercept]"  
##  [29] "r_BIRD[nest-box6,Intercept]"   "r_BIRD[nest-box7,Intercept]"  
##  [31] "r_BIRD[nest-box8,Intercept]"   "r_BIRD[nest-box9,Intercept]"  
##  [33] "r_BIRD[other1,Intercept]"      "r_BIRD[other10,Intercept]"    
##  [35] "r_BIRD[other2,Intercept]"      "r_BIRD[other3,Intercept]"     
##  [37] "r_BIRD[other4,Intercept]"      "r_BIRD[other5,Intercept]"     
##  [39] "r_BIRD[other6,Intercept]"      "r_BIRD[other7,Intercept]"     
##  [41] "r_BIRD[other8,Intercept]"      "r_BIRD[other9,Intercept]"     
##  [43] "r_BIRD[tree1,Intercept]"       "r_BIRD[tree10,Intercept]"     
##  [45] "r_BIRD[tree2,Intercept]"       "r_BIRD[tree3,Intercept]"      
##  [47] "r_BIRD[tree4,Intercept]"       "r_BIRD[tree5,Intercept]"      
##  [49] "r_BIRD[tree6,Intercept]"       "r_BIRD[tree7,Intercept]"      
##  [51] "r_BIRD[tree8,Intercept]"       "r_BIRD[tree9,Intercept]"      
##  [53] "r_BIRD[inside1,MONTHJan]"      "r_BIRD[inside10,MONTHJan]"    
##  [55] "r_BIRD[inside2,MONTHJan]"      "r_BIRD[inside3,MONTHJan]"     
##  [57] "r_BIRD[inside4,MONTHJan]"      "r_BIRD[inside5,MONTHJan]"     
##  [59] "r_BIRD[inside6,MONTHJan]"      "r_BIRD[inside7,MONTHJan]"     
##  [61] "r_BIRD[inside8,MONTHJan]"      "r_BIRD[inside9,MONTHJan]"     
##  [63] "r_BIRD[nest-box1,MONTHJan]"    "r_BIRD[nest-box10,MONTHJan]"  
##  [65] "r_BIRD[nest-box2,MONTHJan]"    "r_BIRD[nest-box3,MONTHJan]"   
##  [67] "r_BIRD[nest-box4,MONTHJan]"    "r_BIRD[nest-box5,MONTHJan]"   
##  [69] "r_BIRD[nest-box6,MONTHJan]"    "r_BIRD[nest-box7,MONTHJan]"   
##  [71] "r_BIRD[nest-box8,MONTHJan]"    "r_BIRD[nest-box9,MONTHJan]"   
##  [73] "r_BIRD[other1,MONTHJan]"       "r_BIRD[other10,MONTHJan]"     
##  [75] "r_BIRD[other2,MONTHJan]"       "r_BIRD[other3,MONTHJan]"      
##  [77] "r_BIRD[other4,MONTHJan]"       "r_BIRD[other5,MONTHJan]"      
##  [79] "r_BIRD[other6,MONTHJan]"       "r_BIRD[other7,MONTHJan]"      
##  [81] "r_BIRD[other8,MONTHJan]"       "r_BIRD[other9,MONTHJan]"      
##  [83] "r_BIRD[tree1,MONTHJan]"        "r_BIRD[tree10,MONTHJan]"      
##  [85] "r_BIRD[tree2,MONTHJan]"        "r_BIRD[tree3,MONTHJan]"       
##  [87] "r_BIRD[tree4,MONTHJan]"        "r_BIRD[tree5,MONTHJan]"       
##  [89] "r_BIRD[tree6,MONTHJan]"        "r_BIRD[tree7,MONTHJan]"       
##  [91] "r_BIRD[tree8,MONTHJan]"        "r_BIRD[tree9,MONTHJan]"       
##  [93] "prior_Intercept"               "prior_b"                      
##  [95] "prior_sigma"                   "prior_sd_BIRD"                
##  [97] "prior_cor_BIRD"                "lp__"                         
##  [99] "accept_stat__"                 "stepsize__"                   
## [101] "treedepth__"                   "n_leapfrog__"                 
## [103] "divergent__"                   "energy__"

pars <- starling.brm3 %>% get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") %>% na.omit()

starling.brm3$fit %>%
  stan_trace(pars = pars)

The chains appear well mixed and very similar

• stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

starling.brm3$fit %>%
  stan_ac(pars = pars)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

starling.brm3$fit %>% stan_rhat()

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

starling.brm3$fit %>% stan_ess()

Ratios all very high.

starling.brm3$fit %>%
  stan_dens(separate_chains = TRUE, pars = pars)

ggmcmc

The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

## starling.ggs <- starling.brm3 %>% ggs(burnin = FALSE, inc_warmup = FALSE)
## starling.ggs %>% ggs_traceplot()

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

## ggs_autocorrelation(starling.ggs)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

## ggs_Rhat(starling.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

## ggs_effective(starling.ggs)

Ratios all very high.

More diagnostics

## ggs_crosscorrelation(starling.ggs)

## ggs_grb(starling.ggs)

Model validation

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

starling.brm3 %>% pp_check(type = "dens_overlay", nsamples = 100)

The model draws appear to be consistent with the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

starling.brm3 %>% pp_check(type = "error_scatter_avg")

This is not really interpretable

• intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

starling.brm3 %>% pp_check(group = "BIRD", type = "intervals")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(starling.brm2)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- starling.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
starling.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = starling$MASS,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(starling.resids, quantreg = TRUE)

Conclusions:

• the simulated residuals do not suggest any issues with the residuals
• there is no evidence of a lack of fit.

Partial effects plots

brms

conditional_effects

starling.brm3 %>%
  conditional_effects("SITUATION:MONTH") %>%
  plot(points = TRUE)

ggpredict

starling.brm3 %>%
  ggpredict(~ SITUATION * MONTH) %>%
  plot(add.data = TRUE)

ggemmeans

starling.brm3 %>%
  ggemmeans(~ SITUATION | MONTH) %>%
  plot(add.data = TRUE)

fitted_draws

Partial.obs <- starling.brm3$data %>%
  mutate(
    Pred = predict(starling.brm3)[, "Estimate"],
    Resid = resid(starling.brm3)[, "Estimate"],
    Obs = Pred + Resid
  )

starling.brm3 %>%
  fitted_draws(newdata = starling, re_formula = NA) %>%
  median_hdci() %>%
  ggplot(aes(x = SITUATION, y = .value, color = MONTH)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper)) +
  geom_line() +
  geom_point(
    data = Partial.obs, aes(y = Obs, x = SITUATION, color = MONTH),
    position = position_nudge(x = 0.1)
  ) +
  geom_point(
    data = starling, aes(y = MASS, x = SITUATION, color = MONTH), alpha = 0.2,
    position = position_nudge(x = 0.05)
  )

Model investigation

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

starling.brm3 %>% summary()

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: MASS ~ MONTH * SITUATION + (MONTH | BIRD) 
##    Data: starling (Number of observations: 80) 
##   Draws: 3 chains, each with iter = 1000; warmup = 200; thin = 5;
##          total post-warmup draws = 480
## 
## Group-Level Effects: 
## ~BIRD (Number of levels: 40) 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)               0.74      0.60     0.02     2.23 1.00     2061
## sd(MONTHJan)                1.34      1.08     0.05     3.84 1.00     1869
## cor(Intercept,MONTHJan)     0.03      0.58    -0.94     0.96 1.00     2314
##                         Tail_ESS
## sd(Intercept)               1992
## sd(MONTHJan)                1968
## cor(Intercept,MONTHJan)     2236
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     79.13      1.14    77.01    81.39 1.00     2354
## MONTHJan                       8.53      1.53     5.48    11.53 1.00     2545
## SITUATIONnestMbox              0.38      1.59    -2.75     3.42 1.00     2577
## SITUATIONother                -3.48      1.59    -6.53    -0.39 1.00     2489
## SITUATIONtree                  4.15      1.60     1.06     7.27 1.00     2399
## MONTHJan:SITUATIONnestMbox     2.02      2.19    -2.24     6.44 1.00     2546
## MONTHJan:SITUATIONother        0.05      2.19    -4.24     4.49 1.00     2389
## MONTHJan:SITUATIONtree        -0.95      2.17    -5.25     3.21 1.00     2474
##                            Tail_ESS
## Intercept                      2335
## MONTHJan                       2327
## SITUATIONnestMbox              2411
## SITUATIONother                 2329
## SITUATIONtree                  2259
## MONTHJan:SITUATIONnestMbox     2173
## MONTHJan:SITUATIONother        2370
## MONTHJan:SITUATIONtree         2428
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     3.97      0.41     3.17     4.81 1.00     2226     2145
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• the estimated Mass of starlings nesting Inside in November is 79.13
• starlings that nest Inside were found to increase their Mass from November to January by 8.53 grams.
• starling that nest in nest boxes are on average 0.38 grams heavier in November than those that nest inside.
• starling that nest in other and trees are on average 3.48 grams lighter and 4.15 heavier respectively in November than those that nest inside.
• there is relatively little variation in Mass between individual birds compared to the variation in Mass within Birds.

tidyMCMC

starling.brm3$fit %>%
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )

Conclusions:

see above

as_draws_df (posteriors)

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

starling.brm3 %>% get_variables()

##   [1] "b_Intercept"                   "b_MONTHJan"                   
##   [3] "b_SITUATIONnestMbox"           "b_SITUATIONother"             
##   [5] "b_SITUATIONtree"               "b_MONTHJan:SITUATIONnestMbox" 
##   [7] "b_MONTHJan:SITUATIONother"     "b_MONTHJan:SITUATIONtree"     
##   [9] "sd_BIRD__Intercept"            "sd_BIRD__MONTHJan"            
##  [11] "cor_BIRD__Intercept__MONTHJan" "sigma"                        
##  [13] "r_BIRD[inside1,Intercept]"     "r_BIRD[inside10,Intercept]"   
##  [15] "r_BIRD[inside2,Intercept]"     "r_BIRD[inside3,Intercept]"    
##  [17] "r_BIRD[inside4,Intercept]"     "r_BIRD[inside5,Intercept]"    
##  [19] "r_BIRD[inside6,Intercept]"     "r_BIRD[inside7,Intercept]"    
##  [21] "r_BIRD[inside8,Intercept]"     "r_BIRD[inside9,Intercept]"    
##  [23] "r_BIRD[nest-box1,Intercept]"   "r_BIRD[nest-box10,Intercept]" 
##  [25] "r_BIRD[nest-box2,Intercept]"   "r_BIRD[nest-box3,Intercept]"  
##  [27] "r_BIRD[nest-box4,Intercept]"   "r_BIRD[nest-box5,Intercept]"  
##  [29] "r_BIRD[nest-box6,Intercept]"   "r_BIRD[nest-box7,Intercept]"  
##  [31] "r_BIRD[nest-box8,Intercept]"   "r_BIRD[nest-box9,Intercept]"  
##  [33] "r_BIRD[other1,Intercept]"      "r_BIRD[other10,Intercept]"    
##  [35] "r_BIRD[other2,Intercept]"      "r_BIRD[other3,Intercept]"     
##  [37] "r_BIRD[other4,Intercept]"      "r_BIRD[other5,Intercept]"     
##  [39] "r_BIRD[other6,Intercept]"      "r_BIRD[other7,Intercept]"     
##  [41] "r_BIRD[other8,Intercept]"      "r_BIRD[other9,Intercept]"     
##  [43] "r_BIRD[tree1,Intercept]"       "r_BIRD[tree10,Intercept]"     
##  [45] "r_BIRD[tree2,Intercept]"       "r_BIRD[tree3,Intercept]"      
##  [47] "r_BIRD[tree4,Intercept]"       "r_BIRD[tree5,Intercept]"      
##  [49] "r_BIRD[tree6,Intercept]"       "r_BIRD[tree7,Intercept]"      
##  [51] "r_BIRD[tree8,Intercept]"       "r_BIRD[tree9,Intercept]"      
##  [53] "r_BIRD[inside1,MONTHJan]"      "r_BIRD[inside10,MONTHJan]"    
##  [55] "r_BIRD[inside2,MONTHJan]"      "r_BIRD[inside3,MONTHJan]"     
##  [57] "r_BIRD[inside4,MONTHJan]"      "r_BIRD[inside5,MONTHJan]"     
##  [59] "r_BIRD[inside6,MONTHJan]"      "r_BIRD[inside7,MONTHJan]"     
##  [61] "r_BIRD[inside8,MONTHJan]"      "r_BIRD[inside9,MONTHJan]"     
##  [63] "r_BIRD[nest-box1,MONTHJan]"    "r_BIRD[nest-box10,MONTHJan]"  
##  [65] "r_BIRD[nest-box2,MONTHJan]"    "r_BIRD[nest-box3,MONTHJan]"   
##  [67] "r_BIRD[nest-box4,MONTHJan]"    "r_BIRD[nest-box5,MONTHJan]"   
##  [69] "r_BIRD[nest-box6,MONTHJan]"    "r_BIRD[nest-box7,MONTHJan]"   
##  [71] "r_BIRD[nest-box8,MONTHJan]"    "r_BIRD[nest-box9,MONTHJan]"   
##  [73] "r_BIRD[other1,MONTHJan]"       "r_BIRD[other10,MONTHJan]"     
##  [75] "r_BIRD[other2,MONTHJan]"       "r_BIRD[other3,MONTHJan]"      
##  [77] "r_BIRD[other4,MONTHJan]"       "r_BIRD[other5,MONTHJan]"      
##  [79] "r_BIRD[other6,MONTHJan]"       "r_BIRD[other7,MONTHJan]"      
##  [81] "r_BIRD[other8,MONTHJan]"       "r_BIRD[other9,MONTHJan]"      
##  [83] "r_BIRD[tree1,MONTHJan]"        "r_BIRD[tree10,MONTHJan]"      
##  [85] "r_BIRD[tree2,MONTHJan]"        "r_BIRD[tree3,MONTHJan]"       
##  [87] "r_BIRD[tree4,MONTHJan]"        "r_BIRD[tree5,MONTHJan]"       
##  [89] "r_BIRD[tree6,MONTHJan]"        "r_BIRD[tree7,MONTHJan]"       
##  [91] "r_BIRD[tree8,MONTHJan]"        "r_BIRD[tree9,MONTHJan]"       
##  [93] "prior_Intercept"               "prior_b"                      
##  [95] "prior_sigma"                   "prior_sd_BIRD"                
##  [97] "prior_cor_BIRD"                "lp__"                         
##  [99] "accept_stat__"                 "stepsize__"                   
## [101] "treedepth__"                   "n_leapfrog__"                 
## [103] "divergent__"                   "energy__"

starling.draw <- starling.brm3 %>%
  gather_draws(`b.Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE)
starling.draw

We can then summarise this

starling.draw %>% median_hdci()

starling.brm3 %>%
  gather_draws(`b_Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE) %>%
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = .variable,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)

starling.brm3 %>%
  gather_draws(`.Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE) %>%
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_halfeye(aes(x = .value, y = .variable)) +
  theme_classic()

bayesplot

starling.brm3$fit %>% plot(type = "intervals")

tidy_draws

This is purely a graphical depiction on the posteriors.

starling.brm3 %>% tidy_draws()

spread_draws

starling.brm3 %>% spread_draws(`.*Intercept.*|b_SITUATION.*|b_MONTH.*`, regex = TRUE)

posterior_samples

starling.brm3 %>%
  posterior_samples() %>%
  as_tibble()

\(R^2\)

starling.brm3 %>%
  bayes_R2(re.form = NA, summary = FALSE) %>%
  median_hdci()

starling.brm3 %>%
  bayes_R2(re.form = ~ (1 | BIRD), summary = FALSE) %>%
  median_hdci()

starling.brm3 %>%
  bayes_R2(re.form = ~ (MONTH | BIRD), summary = FALSE) %>%
  median_hdci()

Further investigations

brms

starling.brm3 %>%
  emmeans(~SITUATION) %>%
  pairs()

##  contrast            estimate lower.HPD upper.HPD
##  inside - (nest-box)    -1.41    -3.783     1.182
##  inside - other          3.48     0.759     5.881
##  inside - tree          -3.68    -6.203    -1.037
##  (nest-box) - other      4.86     2.321     7.406
##  (nest-box) - tree      -2.30    -5.026     0.263
##  other - tree           -7.14    -9.789    -4.419
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

##  contrast  estimate lower.HPD upper.HPD
##  Nov - Jan    -8.81     -10.6     -7.08
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

## SITUATION = inside:
##  contrast  estimate lower.HPD upper.HPD
##  Nov - Jan    -8.56     -11.4     -5.38
## 
## SITUATION = nest-box:
##  contrast  estimate lower.HPD upper.HPD
##  Nov - Jan   -10.54     -14.1     -7.22
## 
## SITUATION = other:
##  contrast  estimate lower.HPD upper.HPD
##  Nov - Jan    -8.55     -11.8     -4.91
## 
## SITUATION = tree:
##  contrast  estimate lower.HPD upper.HPD
##  Nov - Jan    -7.58     -11.0     -4.29
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

Compare specific sets of SITUATION within each MONTH

levels(starling$SITUATION)

## [1] "inside"   "nest-box" "other"    "tree"

## MONTH = Nov:
##  contrast        estimate lower.HPD upper.HPD
##  SITUATION.Comp1    3.666     0.855      6.49
##  SITUATION.Comp2    1.544    -1.131      4.11
## 
## MONTH = Jan:
##  contrast        estimate lower.HPD upper.HPD
##  SITUATION.Comp1    4.610     1.401      7.69
##  SITUATION.Comp2    0.552    -2.590      3.46
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

References